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In this work we present a new method to compute the delays of delay differential 
equations (DDEs), such that the DDE has a purely imaginary eigenvalue. For delay 
differential equations with multiple delays, the critical curves or critical surfaces 
in delay space (that is, the set of delays where the DDE has a purely imaginary 
eigenvalue) are parameterized. We show how the method is related to other works 
in the field by treating the case where the delays are integer multiples of some delay 
value, i.e., commensurate delays. 

, ^ , The parametrization is done by solving a quadratic eigenvalue problem which is 

constructed from the vectorization of a matrix equation and hence typically of large 
size. For commensurate delay differential equations, the corresponding equation is 
a polynomial eigenvalue problem. As a special case of the proposed method, we find 
a closed form for a parameterization of the critical surface for the scalar case. 

We provide several examples with visualizations where the computation is done 
\q with some exploitation of the structure of eigenvalue problems. 

O 
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1 Introduction 



Some phenomena in engineering, physics and biology can be accurately de- 
scribed with delay- differential equations (DDEs), sometimes referred to as 
time-delay systems (TDS). The delay in these models typically stem from 
modelling of phenomena like transmission, transportation and inertia. For in- 
stance, in engineering it is often desirable to design a controller for a system 
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where physical limitations makes the current state of the system unavailable 
for measurement. Instead, the feedback in the controller is done with an old 
state of the system. This is a typical case where a delay-differential equation 
may be used to model the behavior. See [23] for more examples. In this paper 
we consider linear, time- invariant, n-dimensional delay-differential equations 
with m delays, 

( m 

x(t) = A x(t) + J2A k x(t-h k ),t>0 

k=i (1) 

x (t) = <f(t),te [-h m ,o] 

with x : [-hm, oo) ^ W 1 and h k G R+, A k G R nxn for k = 0, . . . , m. 

For most applications it is desirable that the solution is stable independent of 
the initial condition, i.e., that the solution goes asymptotically to zero when 
time goes to infinity. This and many other property can be determined from 
the solutions of the characteristic equation. The characteristic equation of ([!]) 
is 

M(s)v := \-sI n + A + p^A k e- hkS ^j v = Q, \\v\\ = 1, (2) 

where v G C n is called eigenvector and s G C an eigenvalue. The set of all 
eigenvalues is called the spectrum. 

Similar to ordinary differential equations (ODEs), a DDE is exponentially 
stable if and only if all eigenvalues lie in the open left complex half-plane, 
because the only clustering point of the real part is minus infinity. An essential 
difference is that, unlike ODEs, the spectrum of DDEs contains a countably 
infinite number of eigenvalues. Fortunately, it can be proven (c.f. [15]) that 
there are only a finite number of eigenvalues to the right of any vertical line in 
the complex plane, e.g., there are only finite number of unstable eigenvalues. 

It is of particular interest to determine for what choices of the delays h±, . . . , h m 
the DDE ([I]) is stable. This set of hi, . . . , h m is referred to as the stability region 
in delay-parameter space. From continuity it is clear that, at the boundary of 
any stability region in the delay-parameter space, the DDE has at least one 
purely imaginary eigenvalue. These critical delay values (single delay), critical 
curves (two delays) or critical surfaces (more delays) are hence important for 
a complete stability analysis 1 . In this work we focus on the computation of 
critical curves and surfaces. 



In Section 2.3 we present a new method to parameterize such curves and sur- 



faces. An important step of the method is to compute the eigenvalues of a 



1 The critical curves (surfaces) are called offspring curves and kernel curves in [30] , 
Hopf bifurcation curves (surfaces) in [13] and crossing delays (curves) in [22] and 
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large quadratic eigenvalue problem, c.f. [[33], such that the eigenvalue is on the 
unit circle. The matrices of the eigenvalue problem are of size n 2 x n 2 and 
hence large even for DDEs of moderate size. The size of the matrices will 
cause a computational bottleneck and we therefore make some initial notes on 
how an implementation can make use of the structure of the problem of the 
corresponding linearized eigenvalue problem. In Section 2A we provide a new 
interpretation of the results by Chen, Gu and Nett [6], where the commensu- 
rate case is treated. The critical delays for a system with commensurate delays 
can be computed from the eigenvalues of a polynomial eigenvalue problem of 
degree 2m. 



The large number of approaches in the literature yielding delay-dependent 
stability conditions can be classified into two main branches. Most results are 
either based on the construction of a Lyapunov-Krasovskii functional, or (as 
here) based on a discussion of the imaginary eigenvalues. 



Both types of approaches have advantages and disadvantages. For instance, 
the methods in [S], [21] and the methods summarized in [TT] are examples 
where a Lyapunov-Krasovskii functional is used to construct sufficient but 
conservative delay-dependent stability conditions formulated as linear matrix 
inequalities (LMIs). An advantage of these types of approaches is that the 
conditions can be elegantly formulated with LMIs which allow the automatic 
application in engineering software. A disadvantage is that the results can be 
conservative and may impose unnecessarily strong constraints. Moreover, the 
treatment of LMIs of large dimension may be computationally cumbersome. 



The approaches based on imaginary eigenvalues typically yield exact results 
but can not be compactly formulated. Early works (and some more recent 
works) with a discussion of the imaginary eigenvalues were limited to special 
cases. For instance, the trigonometric conditions by Nussbaum [26j Chapter 
III] are for scalar problems with only two delays, the conditions by Cooke and 
Grossman [7] for second order DDEs with a single delay, the geometric char- 
acterization by Hale and Huang in [T3] for first order (scalar) systems with 
two delays, the discussion using the quasi-polynomial by Beretta and Kuang 
in [3] for single-delay DDEs where the coefficients are dependent on the delay, 
the analysis of two-delay DDEs of a special form in [12] (further discussed in 
Example [l]), the analysis of third order DDEs with two delays by Sipahi, et 
al in [3U] based on the earlier methods by Thowsen [32], Rekasius [27] and 
Hertz [TB], the analysis based on a Pontryagin's results for quasipolynomials 
for second order (scalar) DDEs by Cahlon et al in [5], the explicit trigonomet- 
ric analysis of scalar two-delay DDEs by Belair et al in [2]. From this large 
(but incomplete) list of works, we conclude that an analysis based on imag- 
inary eigenvalues is indeed accepted as a natural way to construct stability 
conditions. 
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Recent focus in this field is towards neutral DDEs, distributed DDEs and con- 
ditions which can be efficiently and accurately computed on a computer. The 
importance of the possibility to apply the conditions using a computer should 
not be underestimated. In fact, the conditions from most of the methods men- 
tioned above are implicitly or explicitly formulated using the coefficients in 
the characteristic polynomial with exponential terms (sometimes referred to 
as a quasipolynomial). It is well known from numerical analysis of eigenvalue 
problem, that this (nominal) representation of the characteristic equation is 
not (numerically) stable, i.e., the eigenvalues can be (and are often) very sen- 
sitive with respect to perturbations (e.g. rounding errors) in the coefficients. 
This is important, since it is not possible to represent a polynomial on a com- 
puter using its coefficients without introducing small errors. For this reason, 
the methods mentioned above applied to systems of order (say) 4 or more 
typically do not give numerical results which can be used in practice. This 
motivates the recent works on imaginary eigenvalues using a formulation of 
the matrices in the characteristic equation which is generally believed to scale 
better with the dimension of the DDE. The following works represents the 
class of conditions formulated as generalized eigenvalue problems, Chen, et 
al [6 J for commensurate systems and derivative works [TUPo] . Louisell [TH] for 
single delay neutral systems. In the conditions of these matrix-pencil methods, 
the eigenvalues of a large generalized eigenvalue problem must be computed. 
Similar to the method discussed here, the matrices are constructed using the 
Kronecker product. 

We also mention the related work by Ergenc, et al [E] which is partially based 
on reason with matrices, but also uses the coefficients in the quasipolynomial. 
The method presented here has similarities with this work. 

Finally, we note that there are methods to numerically determine parts of the 
spectra, e.g. using a linear multistep discretization of the solution operator [31] 
or discretization of the infinitesimal generator of the semigroup corresponding 
to the DDE [I]. These methods can be applied for a grid of points in delay- 
space. If the grid is fine enough, a numerical stability condition is computed. 

Apart from the standard works in the field [TTf2"3"] . we refer the reader to the 
references in the introduction of [18] for a well balanced overview of modern 
delay-dependent stability conditions. 



2 Results 



In order to clearly state the results we start by presenting a motivating example 



in Section 2.1 The example describes the problem as well as the idea behind 
the method. 
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The main result consists of the method to parameterize the critical curves 
(surfaces) in Section 2.3 and some further interpretation and relation to other 
publications for the case that the delays are commensurate in Section 2.4 The 
computationally dominating part of both of the methods is the determination 
of unit eigenvalues of quadratic or polynomial eigenvalue problems of large 
dimension. 



2. 1 Motivation 



We now describe the general idea behind the method presented in Section 2.3 
We do this by first considering the scalar two-delay case and describe how 
the derivation must be changed to hold for multiple dimensions and multiple 
delays. 

We wish to find purely imaginary eigenvalues, say s = iu> with corresponding 
eigenvector v, i.e., M.(iu)v = 0. In order to clearly motivate the results we 
discuss the two delay case and generalize and formalize the discussion in the 
following subsections. The characteristic equation for a DDE with two delays 
is 

M(iuj)v = (-iul + A + A 1 e~ ihlUJ + A 2 e~ lh2UJ ) v = 0. 

Except for some special cases, the points (in delay space) of interest, i.e., points 
which causes the DDE to have imaginary eigenvalues, correspond to curves. 
We wish to parameterize these critical curves, and pick (p := h±u as a free 
parameter. By treating ip as a known value we form a condition on h 2 and uj 
(such that we can also compute hi). We will see later that this choice of free 
parameter is valid for most cases. The characteristic equation is 27r-periodic in 
tp, and it will be enough to let ip run the whole span [—it, ir], i.e. for each choice 
of <p in this interval we will find some points on the critical curves and if we 
let ip run the whole interval we will find all critical points. The characteristic 
equation is such that we wish to find all uj G R and z G dD such that 

M(iu)v = (-iul + A + A x e~ iLp + A 2 z) v = 0, v ^ (3) 

where dD denotes the unit circle. If we first consider the scalar case, i.e., A Q = 
a ,Ai = a ly A 2 = a 2 G K the equation corresponds to two scalar conditions 
(say real and imaginary parts) and we can eliminate either uj or z. In the 
approach presented here, we eliminate tu by forming the sum of ^ and its 
complex conjugate, i.e., 

= 2a + a\e~ l{f> + a 2 z + aie 1 ^ + a 2 z = a 2 z + 2ao + 2ai cos(^) + a 2 z. 

Multiplying with z yields the quadratic equation, 

a 2 z 2 + 2z(a + a% cos(^)) + a 2 = 0, 
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since zz = 1. It has the two solutions 

— (do + ai cos(v9)) ± iJa% — (oq + «i cos(<^))' 



a 2 



assuming 02 7^ 0. We can now compute to by inserting z into (|3]) and rearrange 



the terms, i.e., 



iu = a + a^e ltp + a 2 z = i ( —a\ sin(y?) ± y a\ — (oq + Oi cos(y?)) 2 
Since z = e~ lh2UJ and = /iicu, a parametrization of the critical curves is 

V 



h(<p) 



hi 
hi, 



( 



—ai sin((/p)±-y/ a| — (ao+ai cos((/p)) 2 
— Arg z+2<pr 



(4) 



\ -ai sin(v3)±y' a?,-(ao+ai cos(</?)) 2 / 



where Arg 2; = ±sign(a,2)acos ^_ a o+ a i^ os (^) ^ f or ari y p^q £ % and (/? G [— 7r, 7r]. 
Note that the signs must be matched, i.e., for each choice of the free parameter 
(p and branches p and q, the parameterization has two critical delays. 

There are other approaches based on other parameterizations. For instance, 
the analysis of Gu et al [12] can be seen as a parameterization with weRas 
free parameter. They reach the result that if the characteristic equation can 
be rewritten to the form 

l + ai (s)e- h ^ + a 2 (s)e- h * 8 = 0, 

where a\ and a 2 are rational functions, a parameterization of the critical curves 
are given by 



Arg(oi(io;)) + (2p-l)7r±0 1 
hi = , v± = acos 



to 



h 2 = , 2 = acos 



1 + 


ai(iuj)\ 2 - 


a 2 {iu)\ 2 


1 + 


2\ai(iu) 
a 2 (iuj)\ 2 - 


1 

ai{iuj)\ 2 


2 


a 2 (iu)\ 



We also note that the different types of possible critical curves and other 
properties are classified in [14J. 

Example 1 (Different parameterizations) In order to show the differ- 
ence between ^ and the parameterization in JTIf . we construct the parame- 
terizations of the critical curves for the case that qq = a± — —1, a 2 = —\. 
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According to Q, the critical curves are given by, 



( ip+2pir \ 

sin(yj)± | \J 5+8 cos(ip)+4 cos 2 (ip) 
±acos(— 2— 2 cos(i/?))+2(j'7r 
\ sin(</?)± | yj 5+8 cos(i/?)+4 cos 2 (ip) j 



In this example we can find a slightly more elegant parametrization by letting 
x = cos(<^), i.e., 



( 



h(x) 



acos(x)+2p7r 



-sign(a:)N/l-a; 2 ±- v /0.25-(l+a;) 
±acos(-2-2x)+2(j7r 
V -sign(x) v / l-a; 2 ± v /0.25-(l+a;) 2 / 



,x e [-1,1]. 



The critical curves are shown in Figure\l\ The minimum of the 2-norm, i.e. the 
2-norm stability delay-radius, is ||/i(x)||2 ~ 2.896 and taken at x ~ —0.7012 
where u « 1.1139, /ii « 2.1078 and h 2 « 1.9853. 
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Figure 1. Critical curves for Example [T] 

In i/ie context of [12], a\{iio) = j^,a 2 (iuj) = 27i+^j) anc ^ the parameterization 
is given by 



, -ATg(l+iu) + (2 P -l)7T±e i (u) 

hi — 



UJ 



— ( — atan(c<j) ± acos ( -VT+~tJ^H ) + (2p — \)n 

to \ \2 8vl + u; 2 / , 

-Arg (1 + iw) + {2q - 1)tt =f 2 (w) 



— f — atan(cj) =)= acos [ \/l + u 2 j + (2q — 1W 

U V V 4V1 + W 2 / , 
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which represent the same set of curves, but has a surprisingly small obvious 
similarity with the other parameterization Q . 

We now wish to do an analysis similar to the derivation of Q for the general 
case where Aq, A\ and A2 are not scalar. It turns out that a straightforward 
generalization does not yield such explicit results. The reason is that the elim- 
ination of uj must be done in a slightly different way. The first step in the 
discussion above involves taking the sum of ^ and its complex conjugate. 
Since v is typically a complex vector, the sum of ^ and its complex con- 
jugate (|3| will not eliminate uj. Taking the sum of the ^ and its complex 
conjugate transpose is of course also not possible as the dimensions do not fit. 
Instead we take the sum of M.(s)vv* and vv*M.(s)*, where A* denotes the com- 
plex conjugate transpose of A. We believe this is a natural way to generalize 
the elimination of to. Clearly, 

= M(s)w* + vv*M(s)* = 

= A 2 vv*z + {{A + Ate-^vv* + vv*(A% + A\e iLp )) + w*A%z, 

which is a matrix equation. This equation can be rewritten into a quadratic 
eigenvalue problem by vectorization, i.e., stacking the columns of the matrix 
on top of each other. Quadratic eigenvalue problems can be rewritten into a 
generalized eigenvalue problem using a so called companion linearization (c.f 
PTj). As we will see later, if n is not large, there are numerical methods to 
find all z. We can then compute the corresponding uj and the critical delays 
similar to the way done for the scalar case. 

Even though the motivation above is simple, it is not clear that all steps in- 
volved are reversible, i.e., there is no convincing argument to ensure that we 
get all critical delays. For that reason, we formalize the discussion in the sec- 
tions that follow. At the same time we generalize the method to an arbitrary 
number of delays. The generalization is such that if we have m delays we 
have m — 1 free parameters, and need to solve a quadratic eigenvalue problem 
for each choice of the parameters. In order to compare the method with re- 
lated approaches we also discuss the case where the delays are (fixed) integer 
multiples of the first delay, i.e., commensurate delays. 

2.2 Main results 

Our main goal is to formalize and generalize the discussion in the previous sec- 
tion, and to see how the constructed method relates to approaches in the liter- 
ature. In particular, we wish to provide further understanding to the method 
in [B] where a matrix-pencil method is presented for the commensurate case, 
i.e., the case where the delays are integer multiples of the first delay. We do 
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this by first, in a general setting, introducing an operator L. For this opera- 
tor we can formulate an equivalence theorem with the characteristic equation. 
One reason for formulating such a theorem is that the results in [5] are easily 
interpreted from the conditions in the equivalence theorem. Another reason 
is that from the conditions we can construct a new method (similar to the 
motivation) to compute the critical delays for the case that the delays are not 
necessarily commensurate. Hence, both results can be interpreted in terms of 
the conditions in the equivalence theorem. 

Theorem 2 Given s 6 C and v G C n , v*v = 1 the following statements are 
equivalent. 

M(s)v = 

L(vv*,s) = A v*M(s)v = 
where 

L(X, s) := (M(s)) X + X (M(s)*) = 

m 

= (A k Xe~ hkS + XA T k e~ hkS ) - 2XRe s, 

k=0 

and h = for notational convenience. 

Proof. The forward implication is trivial from definitions, i.e., (|2]) and (|7j). 
The backward implication, i.e., (|6])=^(|5]), is clear from the following equality. 

M(s)v = L(vv*, s)v - vv*M(s)*v 

□ 

The reason why we construct L in this way, is that similar to the motivation, 
for the critical case, i.e., s = iuo e iR, s = iu only appears in the exponential 
terms of the expression. 

In the two following sections we rephrase the equivalence theorem such that 
we get a method to compute the critical delays. Because of the equivalence 
theorem the formulation will not miss any critical delays. 

2.3 Free delays 

In the motivation we saw that for two delays we needed to parameterize a 
curve in delay-space. In general, if we have m delays we need to parameterize 
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(5) 
(6) 

(7) 



the (m — 1) dimensional (hyper-)surface. We do that by letting z = e~ lhmUJ 
and represent the other exponential terms by free variables, i.e., tp k = h k ui, 
k = 1, . . . , m — 1. This allows us to rewrite the main expression of Theorem [2] 
to a quadratic eigenvalue problem which can be solved by transforming it to 
an eigenvalue problem using a companion linearization. 

For the scalar case, the quadratic eigenvalue problem can be solved explicitly, 
and we can form an explicit parametrization of the critical surfaces. We also 
verify the results by comparing it to classical results for the single delay, scalar 
delay-differential equations. 

The following theorem follows from the substitution z = e~ lhmUJ ,ip k = h k cu, 
k = 1 , . . . , m — 1 and Theorem [2j 

Theorem 3 The critical delays are given by h m = ~ Axs z + 2 P™ n ; h k = yfc+ j Pfc7r , 
k = 1, . . . , m — 1, where p k G Z, k = 1 . . . m, ip k G [— vr, tt], k = 1, . . . , m — 1, 
z G 3D, uj G M. such that 

z 2 A m vv* + z (j2 A k vv*e- iVk + vv*A T k e i ' Pk \ + vv* A T m = 0, (8) 

\fc=0 / 

Cm-l \ 
A m z + Y, Me~ iVk v, (9) 
fc=0 / 

where ifo = for notational convenience. 

Apart from singular cases (which will be discussed later), this theorem gener- 
ates a parameterization of the critical curves, with <pi, . . . , ip m -i as free vari- 
ables. It is clear that if we can find solutions z, v of ([8]), u and critical delays 
can be computed for a specific choice of (p^, k = 1, ... ,m — 1. Hence, the 
difficulty of applying this theorem is solving (|8]). We solve this by vectorizing 
and rewriting the quadratic equation to a generalized eigenvalue problem. 



The vectorized version of (|8j), i.e., the same equation but the columns of the 
matrix condition stacked on each other, is the equation 



m— 1 



z 2 I® A m + z Y L k {ipk) + A m ®I \u = 0, (10) 



fe=0 



where u = v ® v is the vectorization of vv* and L k ((p k ) '■= I <%> A k e tlfik + A k ® 
Ie lipk , and ® the Kronecker product with the usual meaning. Problems of the 



type (Mz 2 + Gz + K)v = 0, e.g., (10), are known in the literature as quadratic 
eigenvalue problems (c.f. [S3])- They are typically solved by transformation 
into first order form, normally using a so called companion linearization, e.g., 
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here 



/ 

^m— 1 





(11) 



\I® A m E™= L k (ip k 
We can hence find the solutions of (fsb by finding the eigenvalues corresponding 



to (11 ), to which standard methods for generalized eigenvalue problem can be 
applied. There are also other numerical methods to solve quadratic eigenvalue 
problems without linearizing it, e.g., Jacobi-Davidson [3T] and second order 



Arnoldi {!]. In this paper we focus on the companion linearization (11). 



With these methods it possible (unless n is large) to find the eigenpairs of 



(11) to sufficient accuracy on a computer. Because it is not possible (in gen- 
eral) to solve eigenvalue problems without errors, we will not get exactly a 
vector corresponding to a rank-one matrix. However, if the eigenvalue is iso- 
lated we can find an approximation close to the rank-one matrix. We pick the 
approximation by choosing the rank-one matrix corresponding to the singular 
value which is not close to zero. All but one singular value should be close to 
zero, because the numerical rank is one. The accuracy of this choice of the 
approximation can be characterized by the following simple result. 

Lemma 4 Let (z, x) be an eigenpair of the quadratic eigenvalue problem (z 2 M+ 
zG + K)x = and x = v <g> v, \\x\\ = 1, \z\ — 1. Say x=x+y=u®u+q is 
an approximation of x, then the sine of the angle between the approximation 
u and the vector v is bounded by 



sin(u,«)| < v/||?/||(||M + /|| + ||G|| + ||iq) + ||g||. 



Proof. Multiplying the characteristic equation from the left with x* and 
adding to both sides, yields 

x*xz 2 = x*(z 2 (M + I) + zG + K)x. 

We exploit that x*x = (u (2) u + q)*(v (2) v) = u*vu*v + q*x from which we 
deduce that 

z 2 \u*v\ 2 = z 2 (u*vu*v) = x*{z 2 (M + T) + zG + K)x - z 2 q*x = 

= x*(z 2 (M + I) + zG + K)x + y*(z 2 (M + I) + zG + K)x - z 2 q*x = 

= z 2 + y*(z 2 (M + I) + zG + K)x - z 2 q*x, 

i.e., 

\u*v\ 2 = l + y*((M + I) + z~ x G + Kz~ 2 )x - q*x 

and 

| u * T ,|2>i_|| y ||(||M + /|| + ||G?|| + ||li'||)-|M|, 
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where we used that \z\ = 1. Finally, the result follows from the fact that 
1 — \u*v I 2 . □ 



In rough terms, the result is that if the error of the eigenvector, denoted by y, 
is small and the distance to the Hermitian rank-one matrix q is small, then the 
angle between v and its' approximation u is small. The lemma also indicates 
the (unfortunate) fact that if the error of the eigenvector is e mac h the accuracy 
of the approximation u is essentially y/e ma ch- 

If we have multiple eigenvalues, the invariant subspace must be sought for 
rank-one matrices. This typically unusual case will not be treated in detail. 

Moreover, we are only interested in eigenvalues z on the unit circle. Because of 
rounding errors, we classify an eigenvalue as a unit-eigenvalue if its modulus 
differs from 1 less than some tolerance, say ^e mac h- 



It is also worth noting that (10) is a quadratic eigenvalue problem and can 
have (so called) infinite eigenvalues. For instance, if A m = the equation is 
independent of z G C, i.e., if the linear term is zero it is fulfilled for any z G dD. 
The interpretation of the theorem is then that the criticality is independent 
of h m , and one of the free parameters is no longer free. Since this can only 
occur if A m is singular, the problems of infinite eigenvalues can be avoided by 
reordering of the matrices (assuming at least one matrix is nonsingular). We 
characterize the phenomena of infinite eigenvalues with the following example. 



Example 5 (Singular case) Consider the DDE with 

( 




A) = , A 1 





With this example we wish to show how an infinite eigenvalue of (10) can be 
interpreted, and how this appears in Theorem [3| For this particular example 
it is actually possible to circumvent the problem with the infinite eigenvalue 
simply by switching Ax and A 2 . Here, we will not approach the problem that 
way, because that type of reordering of the matrices is not always possible. 



If (10) has an infinite eigenvalue, it must have a corresponding eigenvector 
such that the quadratic term disappears, i.e., here u = (0, 1) T £g> (0, 1) T = 
(0, 0, 0, 1) T is the only eigenvector of the correct form. For this v, the quadratic 
and the constant term in J8| are zero, and ^ is independent of z. Since z ^ 0, 
v is a valid solution of (pF only if 

A lV v*e~ iipl + vv*A^e itpi = 0. (12) 
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This means that the only choice of ipi which generates a critical delay (corre- 
sponding to this v) is when (12) is fulfilled. Here (12) and hence Q8J) is fulfilled 
only if e = and for ip\ = — | + 2pir with the corresponding critical curve 
hi — — f + 2p7r (for any h 2 ). The critical curves are plotted for three choices 
of e in Figure^ From the figure it is easy to identify that for e = there are a 
vertical critical curves, i.e., curves which are independent of h 2 corresponding 
to the contribution of the infinite eigenvalue. If e > there are no vertical 
lines, no infinite eigenvalues and the critical curves are characterized by the 
finite eigenvalues. 




Figure 2. Critical curves for Example [5] 



For the scalar case, the quadratic eigenvalue problem reduces to a quadratic 
equation. The theorem can then be simplified to 

Corollary 6 For the case that ij. = 6 K, k — 0, . . . , m, the critical delays 
are given by 



p k + 2p k 7t 
hk = , k — 1 m — 1 



UJ 



Tsign(a m )acos ( ^ fc=0 °^ c ° s(yfc) ) + 2p m n 
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u = -^2a k sin(y3 fc ) ± 

k=0 



CLk COs(ipk) 



\k=0 



where p k G Z, k = 1 . . . m, ip k G [— n, ir], k = 1, . . . , m — 1. 



13 



It is worth noting that the formulas in the literature for the scalar single delay 
case normally assume ai < — [ ct Q | , because otherwise the DDE is (un) stable 
independent of delay or the corresponding delay-free DDE is not stable. From 
Corollary [6j the critical delays for the single delay scalar DDE are given by 

— sign(ai)acos(— — ) + 2p7t 

h = , p , 

which reduces to classical results if a\ < 0, c.f. [7] or J23J Section 3.4.1] 



2.4 Commensurate delays 



For the case that the delays hk are integer multiples of the smallest delay, 
say h, Theorem [2] can be restated such that it is very similar to the re- 
sults by Chen, Gu and Nett in [6]. For hk = n^h, the problem is to find 
the critical delay h such that the choice of the delays in the direction h = 
(hi, hi, ■ ■ ■ , h m ) = (riih, n 2 h, . . . , n m h), where ni, . . . ,n m G N, generates a 
purely imaginary eigenvalue. In other words, we fix a "rational" direction in 
delay-space with (n%, . . . , n m ) and search for critical delays along this direc- 
tion. 

-ihu) 



The following theorem follows from Theorem |2| by letting z = e 

Theorem 7 The critical delays h e IR + corresponding to the direction in 
delay space defined by 

h = (hi, ...,h m ) = (hm, hn 2 , hn m ), 

are given by 

— Arg z + 2pn 

h = , 

for any p G Z and for v G C n , v*v = 1, u G R, Z G dD fulfilling 

m 

£ (A k vv*z n " + vv*A T k z~ nk ) = 0, (13) 

k=0 

and 

where no = for notational convenience. 

Analogously to the previous section, if it is possible to find the solutions z 
and v of (13) the other expressions explicitly yield uj and h. Without loss of 
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generality we let n m = maxke[i,...,m] n k- After vectorizing the matrix equation 
(13) we find that 



fe=0 



n m +n k 



+ A k ®Iz 



n m -n k 



where aeC" is the vectorization of vv*, i.e., u 
the form 



N 



u 



0. 



u = (14) 
v ®v. This equation is of 

(15) 



k=0 



which in the literature is known as a polynomial eigenvalue problem. Similar 
to quadratic eigenvalue problems, the most common way to solve polyno- 
mial eigenvalue problems is by companion linearization, which is analyzed 



and generalized in [21] and |2QJ. For instance, the eigenvalues of (15) are the 
eigenvalues corresponding to the generalized eigenvalue problem 



(i 



\ 



B 



( 



w 



nJ 







-Bn B 



N-2 



\ 

I 

■Bn-i ) 



w. 



(16) 



where w = (u T , zu T , z 2 u T , ■ ■ ■ z N 1 u T ) T . By selecting B k according to the co- 
efficients in (14), we can compute all solutions z,v to (13) by solving the 



eigenvalue problem (16). Other methods, such as the Jacobi-Davidson method 



could be applied directly to the polynomial eigenvalue problem without lin- 
earizing it, c.f., [3T] . 

The Hermitian rank-one matrix can be chosen similar to the choice in the 
previous section, i.e., the principal vector in the singular value decomposition, 
since the accuracy result in Lemma [4] generalizes to the polynomial case. 

Lemma 8 Let (z,x) be an eigenpair of the polynomial eigenvalue problem 

1. Say x = x + y = u®u + qisan 



(15) and x 



v, \\x\ 



1, z 



approximation of x, then the sine of the angle between the approximation u 
and the vector v is bounded by 



sin(w, v)\ < » 



Ml fi + EllW) + N 

V k=0 ) 



Proof. The proof is analogous to the proof of Lemma |4j 



□ 



Remark 9 For the case h k = hk, the companion form (16) is very similar 
to the eigenvalue problem occurring in f^. However, in that context it is not 
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recognized that the eigenproblem to be solved is a polynomial eigenproblem and 
that the eigenvector u is the vectorization of an Hermitian rank one matrix. 



2.5 Notes on computation 



We now clarify how we can use Theorem [3] to generate the critical surfaces, 
discuss computational bottlenecks and suggest possible improvements and al- 
ternatives. 

Theorem [3] is an equivalence theorem between cpi, . . . , y? m _i and hi, ... , h m . 
Hence, we can see it as a parameterization of the critical surface. In other 
words, if we let (/?!,... , (p m -i run the whole interval, we will generate all critical 
points. In practice, we typically let the free parameter run over finite number 
of grid points with a grid size small enough such that we can convince ourselves 
of the continuity of the critical curves. 

With this in mind, we state in pseudo-code how to generate the critical curves 
for the two-delay system. 

1. FOR (p = -7T : A : vr 



2. Find eigenpairs (zk,Uk) of (11) 

3. FOR k = 1 : length(z) 

4. IF Zk is on unit circle 

5. Compute Vk such that Uk = vec VkV* k 

6. Compute uok = -iv* k (A 2 z k + Aq + Aie~ lLp ) Vk 

7. Accept critical points (hi,h2) 

_ (p + 2pn 

— j P — Pmax t ■ ■ ■ iPmax 

h _ -Arg z k + 2qn 

) Q — Pmax ; • • • j Pmax 



8. END 

9. END 
10. END 

In step 1, A is the stepsize of the parameter <p. In step 7, p ma x is the number 
of branches which should be included in the computation. Step 7 is not com- 
putationally demanding. We can select p ma x so large that the computation 
contains all relevant branches, say such that all delays smaller than some de- 
lay tolerance are found. This is possible because the delays are monotonically 
increasing or decreasing in the branch parameter. The generalization to more 
than two delays is straighforward. It involves a nesting of the outer iteration 
(step 1) with for-loops of the new free variables ipk and computing the other 
delays in step 7 similar to hi. 
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For DDEs of large or moderate dimension, the main computational effort in 
an implementation of the method is the solving of eigenvalue problems. This 
can be seen from the following remarks on the computational complexity. 

The current general full eigenproblem solvers (for instance the QR-algorithm) 
require, roughly speaking, a computational effort proportional to N 3 to com- 
pute the eigenvalues of a iV x iV-matrix to sufficient accuracy. For Theo- 
rem [3j the companion matrices are of dimension N = 2n 2 and for Theorem [7j 
iV = 2n m n 2 . For both cases the computational effort is hence proportional 
to n 6 , which is a lot for moderate or large system-dimensions n. For sparse 
eigensolvers the complexity is reduced, but the dimension of the system still 
causes a computational bottleneck. 

Remark 10 We note that the very high computational cost for large sys- 
tems, even suggests that instead of parameterizing using m — 1 free variables 
over ip G [— 71", 7r] m ~ 1 , where each evaluation costs n 6 , it may be more effi- 
cient to search ip G [— 7r,7r] m and as an additional constraint check if (or 
how far from) an imaginary eigenvalue is contained in the spectrum of A = 
^0 + Efeli A k e~ tlfk . Computing the spectrum of A has the computational cost 
n 3 , which (for large n) compensates for the additional degree of freedom which 
has to be scanned. We note that some form of iterative method is likely to be 
necessary in order to impose the additional constraint (the imaginary eigen- 
value) in a reliable way. 

This remark was brought to our attention by Michiel Hochstenbach. 

In order to make larger problems tractable we make some initial remarks on 



how to adapt numerical methods for (11). The eigenvalues we wish to find 



have the following special properties we can make use of. 

(1) The eigenvalues of interest lie on the unit circle. 

(2) The matrices resulting from the companion form are sparse. 

(3) The matrices have a very special structure originating from the Kronecker 
construction. 

(4) Only eigenvectors of the form u = vec vv*, i.e. a vectorization of an 
hermitian rank one matrix, are of interest. 

(5) The eigenvalues move continuously with respect to the parameters. 

In this work we only exploit the first two items. The unit eigenvalue property 
is exploited by transforming the unit eigenvalues to the real line, i.e. we make 
the Cayley-transformation z = jz^- That is, if z is a unit eigenvalue of the 
eigenproblem Av = zBv then a is a real eigenvalue of the eigenproblem (A — 
B)v = a[iA + iB)v. The reason for this is that it is believed that finding real 
eigenvalues is an easier problem than finding eigenvalues on the unit circle. 
Here, we apply a method doing rational Krylov scans (c.f. [28], [22]) along the 
real line implemented in the command sptarn in the Matlab control toolbox. 
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However, we note that further exploitation of the properties are necessary if 
we wish to treat large dimensional DDEs. An example of such a structure- 
exploitation is given in [17] . 



3 Examples 



Example 11 (n = 35, m = 3) In order to determine how well the method 



scales with dimension we apply the method in Section 2.5 to the discretization 
of a partial delay- differential equation. We consider the three delay partial 
differential equation, 



' u' t {x, t) = u" x (x, t) + (3(1 + sin(37rx))-u(x, t) 

— k 5(x — xo)u(0, t — hi) 

— K\5{x — xi)u(xi, t — h 2 ) 

— K 2 5(x — x 2 )u(l, t — h s ) 
u x (0,t)=0 

u x (l,t) = 0, 



(17) 



where we pick k = k 2 = 4, K\ = 10, x = 1/3, x\ = 1/2, x 2 = 3/4 and [3 = 10. 



The physical interpretation of equation (17) is the heat equation on a rod with 



length one with heat production over the whole rod causing instability and three 
delayed stabilizing pointwise feedbacks. This system is discretized with central 
difference in space with n equally distributed intervals, yielding a system of the 
formx(t) = A x(t) + A 1 x(t-h 1 ) + A 2 x( y t-h 2 )+A 3 x( y t-h 3 ), where x(t) Gl", 
i.e., a three delay DDE. For this example we pick a step-length such that the 
discretized system is of dimension n = 35. 

The critical surface closest to the origin of the discretized system is plotted in 
Fig. [3| The algorithm computes points on the critical surface, but for visual- 
ization purposes, we connect the points to form a surface. 



The application of the algorithm in Section IL5 with the described exploitation 
for 390 different combinations of the parameters requires 62 minutes on a 
computer running Linux and Matlab 1 on a 2. 4 GHz Intel Pentium 4 processor 
with 512 Mb RAM. 



4 Conclusions 



We study the conditions on the delays for multiple delay time delay systems 
such that the system has an imaginary eigenvalue. This is done by introducing 
a condition using a function L, for which the imaginary eigenvalue condition 
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Figure 3. Boundary of stability region for Example 11 



has a particularly easy form. For the commensurate case the condition can be 
reduced to the computation of the eigenvalues of a matrix similar to that of 
Chen et al [6]. For the general (free) case we propose a new method, which 
involves solving a quadratic eigenvalue problem. 

For the scalar case, the quadratic eigenvalue problem can be explicitly solved 
and we can find a closed explicit expression for the condition on the delays. 
For systems of large dimension we make initial remarks on what properties of 
the problem an adapted numerical scheme should exploit. 

In the examples section we show the efficiency of the method by applying it to 
previously solved examples as well a previously unsolved three-delay example 
of larger dimension. 
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